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Abstract. We revisit an algorithm by Skeel et al. [5, 16] for computing the 

^^ modified, or shadow, energy associated with the symplectic discretization of 

• Hamiltonian systems. By rephrasing the algorithm as a Richardson extrapo- 

,S^ lation scheme arbitrary high order of accuracy is obtained, and provided error 

LJ estimates show that it does capture the theoretical exponentially small drift 

a associated with such discretizations. Several numerical examples illustrate the 

theory. 



^-H 1- Introduction. Numerical simulation of conservative differential equations re- 

^ quires special care in order to avoid introducing non-conservative, or non-physical 

^^ truncation error effects. For Hamiltonian ODEs or Euler-Lagrange equations orig- 

^_ inating from variational principles there exists much evidence [6, 9, 11, 14] that the 

(^—^ proper discretization scheme should be symplectic [7, 9, 11, 18]. In the Hamilton- 

ian case this can be achieved by imposing special conditions on classical methods 
or by methods based on generating functions [7]. In the variational formulation 
^^ symplecticity is achieved by discretizing the action integral and carrying out a dis- 

T-H Crete variation [10]. In some cases these formulations and methods turn out to be 

L| equivalent by the Legendre transformation [8, 10]. 

Focusing on the Hamiltonian side, symplecticity implies that the trajectory 
produced by the numerical algorithm is the exact solution [12] of another, non- 
autonomous "modified" Hamiltonian system close to the original one. Various sta- 
bility results for Hamiltonian ODEs then apply, leading to an understanding of the 
dynamics of such discretizations schemes [6, 9, 11, 15]. Early results on modified 
equations focused on the autonomous part [2, 4, 6, 13, 14] and established that its 
fiow is exponentially close to the numerical trajectory. This work was motivated 
by the bounded error in energy observed in simulations with symplectic schemes. 
The early results are contained in the newer results since the time-dependent part 
is exponentially small due to analyticity. Despite its smallness the non-autonomous 
term excites instabilities through resonances, one consequence being a drift in the 
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modified energy. In simulations requiring millions of steps such as in molecular 
dynamics [9, 17] and celestial mechanics [19] these effects become significant and it 
becomes important to understand and control them. 

Constructing the modified Hamiltonian is equivalent to evaluating many terms 
in the Baker-Campbell-Hausdorff formula, or its continuous analogue [13], a com- 
binatorially complicated task possible only for small systems and to a low order of 
accuracy. Recently, Skeel and coworkers [5, 16] devised a method for numerically 
computing the value of the modified Hamiltonian along the numerical trajectory, 
thus allowing us to track the possible drift in the modified energy. In this paper 
we simplify this method, possibly at the cost of extra storage, and provide expo- 
nentially small error bounds when it is applied in the asymptotic regime. It is then 
used to verify, and justify the theory of modified equation on several test equations 
and methods. 

2. Modified equations. As alluded to in the introduction, the numerical solution 
of an ODE y' — f{y) is interpolated by the exact solution of a modified ODE 
y' — f{%j,t). The modified equation is non-autonomous, but the non-autonomous 
part is exponentially small in the step size. More precisely, given an analytic vector 
field / and a one-step method defined by an analytic mapping '^uji there exists an 
analytic vector field f{y,t), /i-periodic and analytic in t, whose exact flow exactly 
interpolates the numerical trajectory {a;„}, x„_|_i — ^;ij(x„). The construction 
in [12] starts by constructing a modified vector field f{y,t) which is only C°° in 
t whose flow interpolates {x„}. This vector fleld is then transformed by a time- 
dependent coordinate transformation into a vector field / analytic in t. 

The domain of analyticity of / plays an important role in the analysis, and we 
have found it useful to assume that / is analytic for all y in a domain of the form 

Vy := [J{z e C : \y{t) - z\^ < f J = |J{z G C : \^{y{t) - z)\^ < r~J 

t>Q t>0 

for some fy > 0, where y{t) = (t>tj{yo) is the trajectory of the smooth modified 
vector field /. This domain is typically smaller than the domain of analyticity of /, 
and depends on the numerical method. In the following we will use the sup-norm 
\\f\\v ~ sup^gu |/('2^)|oo- With these definitions the main result of [12] in the limit 
/i — > can be formulated as 

Theorem 1. Let "^hj be a one-step method applied to the analytic vector field f, 
and yn+i = ^hjiVn) be the approximations obtained by iterating ^> . Then there 
exists a modified vector field f{y, t) = f(y) + ri (y) + r2{y, t) which is h-periodic in t 
and analytic in {y,t) G V such that its exact flow satisfies y{nh) — ^nfijivo) = ^n- 
In the limit h ^ we have the estimates 



V 



ll/b'<Y^II/l 

F2||t)' ^O — — exp -?7- 



for Q < rj < 1, < 5 < fy. The domain of analyticity of the modified vector field is 

P' = ((z,t) e C^ X C : |3(z - y{t))\oo < Py - S, |5(t - t)\ < ^^ 
I WfWve 

and the norm \\ ■ ||x) is defined by ||/||p' = sup^^ .^^^^p' |/(2,t)|oo- 
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For a Hamiltonian vector field / and a symplectic numerical method [6, 7, 14], 
the modified vector field / is also Hamiltonian [4, 7, 9], with Hamiltonian H = 
H{y) + Gi{y) + G2{y,t) where H, Gi and G2 are the Hamiltonians corresponding 
to the vector fields /, ri and r2, respectively. The change in the modified energy 
along the numerical trajectory therefore satisfies 

where {F,G} := J2T=i ^ij^Pi ~ ^Pj^ij ^^ ^^'^ Poisson bracket. By Theorem 1 this 
drift is very small for small h, thus motivating symplectic methods. 

3. The method of Skeel et al. for computing the modified energy. Skeel 
and coworkers [5, 16] found an ingenious way of evaluating the modified energy H 
at the points {x„} for discretizations based on splittings [3]. Suppose we have an 
Hamiltonian given by H = ^p^M~^p + U{q). An explicit splitting algorithm with 
step size h is given by 

for n = 0,1,2,... 

PO=Pn, qo^Qn 

for s = 1 : 5 



Ps ^Ps-l ~ hasUqiis-i) 
(is = gs-i + hbsM^^ps 



(1) 



end 

Pn+i =Ps, qn+i = qs 

end 

leading to approximations Pn+i = Ps, In+i = Qs when qo = Qm Po = Pn at i„ = nh. 
By choosing the coefficients as , bs appropriately, a method of arbitrary high order 
can be found. The modified Hamiltonian can be found by representing the inner 
loop of (1) as a concatenation of exponential operators [7] 

^h,f{p,q) = e-x.p{-haiUqdp){p, q) exp{hbiM'^pdq) ■ ■ ■ 

eyip{-'hasUqdp) exp{hbsM~^pdq), 

whereby the Baker-Campbell-HausdorfF (BCH) formula is used to find an expres- 
sion so that "^hjix) — exp{hfd){x). 

The approach of Skeel et al. for computing values of the modified energy is to 
append one scalar equation to the numerical integrator, 

for n = 0,1,2,... 

PO=Pn, qQ=qn, /3o=/3n 

for s = 1 : 5 

Ps =Ps-i -hasUq{qs-i) 

Ps = A-l - has[q^^^Uq{qs^i) + 2C/(g,_i)) (2) 

qs ^ qs-i+ hbsM^^ps 
end 

Pn+i = Ps- qn+1 = qs, /3«+l = i^S 
end 
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where (3q = 0. To understand how the modified energy can be recovered from 
{PniQn, Pn}^ notc that by Theorem 1 the numerical trajectory (pn,9n) is exactly 
interpolated by the flow of a Hamiltonian H. The discretization (2) is the dis- 
cretization of Ha — a^H{a^^p,a^^q) where /3 is conjugate to a, the so-called 
homogeneous extension of H{p,q). The discovery in [16] rests on the fact that 
homogeneous extension is a Lie algebra homeomorphism. Thus, since H is con- 
structed by Poisson brackets as in the BCH formula, the modified Hamiltonian for 
(2), Ha, is the homogeneous extension of if, hence the trajectory generated by (2) 
is interpolated by 

q = Hp{p,q,t), 

P' = -Hq{p,q,t), 

^' = q'^H^ip, g, t) + p^Hpip, g, t) - 2H{p, q, t) , 

from which 

H^liq^H, + p^Hp-l3') = li-fp'+p^-q'-n (3) 

The equation for a is removed from (3) since Ha does not depend on /3, the conju- 
gate variable of a, and hence a' = which is solved exactly by the methods we are 
considering. 

Thus the value of H can be computed by finding the derivatives of the interpo- 
lating trajectory (which are not known since we do not have Hp,Hq). In [5, 16] 
estimates of the derivatives are computed using backward difference formulas and 
interpolating polynomials with stored values of {pn,Qn, Pn}- These polynomials 
can be precomputed, but unfortunately the required expressions are very large, and 
they only provide expressions up to order 24. Their method does however have an 
advantage in requiring less stored values than one based on centered differences, 
which might be important if the modified energy is part of the simulation [5] . 

4. Richardson extrapolation. Our suggestion is to use Richardson extrapolation 
in order to avoid the large expressions that arise in the method described in the 
previous section. 

First consider the use of Richardson extrapolation to find the derivative of a 
function, say y, at zero given the function values on a grid. We define the central 
difference approximations 

yjjh) - y{ ~jh) 

2jh 
and compute the Richardson table entries 

Tj^k — Tj-i,k 

We then have by standard results that Tj_j = 2/'(0) + 0{h'^^). In fact, it is straight- 
forward to prove by induction that the Tj^ satisfy 

^ _^^ 2(-ir+Mj)i ^ 

'' " ^ (* - 1)! (fc - 0! (2j - fc + ^)k {j-k + i) '^-'=+^'^ 
where the Pochhammer symbol denotes the falling factorial: 

{n)k = n{n - X){n - 2) . . . (n - fc + 1) = — . 

fc! 



rr, yyj'- y\ j'" ■ , 

^jM = ^71 ' J = 1' 



T,^k+i^T,^k+ ,^' \_,Ji, '\ , fc-l,...,j-l. 
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It follows that the diagonal entries in the Richardson table are given by T^^ = 
DmHiO) where Dm{y) denotes the central-difference approximation to the deriva- 
tive y'(0) using 2m points, defined by 

^ jh[m — j)\[m + J)'. ' 

This approximation satisfies Dmy{0) — y'{0) + 0{h?"^). By choosing the index m 
appropriately, it is possible to find an exponentially accurate approximation for the 
derivative of analytic functions, as stated in the following lemma. 

Lemma 2. Let y{t) be analytic in {i G C : |S(i)| < p}, then there exists an m* 
(which depends on h) and a constant Ci > such that 

|y'(0)-A„.y(0)| <Ci^^^^5^hH||y||^ 

where \\y\\p == sup|cj(^)|<p |y(r)|oo. 

The proof of this result and other results are found in the Appendix. 

5. Computing the modified energy using Richardson extrapolation. Re- 
turning to the computation of the modified energy (3), the derivatives in this 
formula can be computed with Richardson extrapolation using the stored values 
of {pn,(lm Pn\- To Compute the modified energy at t = nh, we define the central 
difference approximation 

^(_-T( ,-. Piin + i)h)-p{(n-j)h) 
2\^ "• ' 2jh 



Tr,i-^[-^'^inh) 



-T, u,q{{n + m-q{{n-m Pi{n + j)h)~(3{{n-j)h) 



p' (nh) 



2jh 2jh 



^ ( t'PiI'+J Pn—j T ^^+J Qn—j Pn-\-j Hn- 



"9n ?r7I + Pn 



2 V " 2jh ^" 2jh 2jh 

for j — 1, . . . and then compute the Richardson table entries as before; 

rpn _ rpn , j,k J-l,fc 1,-1 T _ 1 (c,\ 

The expression T" ■ is a convenient way of computing the approximations and in 
addition it gives a way of estimating the error in the approximation |iJ— Tj_ij_i| « 
\Tjj — Tj-ij-i\^ which is useful for finding a stopping criterion for the extrapolation 
process. 

We mention in passing that accurate values of H might be obtained using Fourier 
series as well [1] , however such methods seem most useful for quasi-periodic motions 
or scattering problems, while the approach taken here seems suitable for a broader 
range of problems. 

The following corollary follows from Lemma 2. 

Corollary 3. Let Pn,qm Pn be given by the numerical scheme then there exists an 
m* such that 

\H{Pn,qn.t^)-T:^,,raA<^ ^''''^l^^\ \qn\imp+\Pn\imp+\mp). 

where p is such that the interpolating trajectory {p{t),q{t), P{t)) is analytic for 

mt)\<p. 
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In Corollary 3 the parameter p related to the domain of analyticity of y{t) is 
undetermined. The following existence lemma will be useful for bounding p. 

Lemma 4 (Domain of analyticity of the solution). Let g(y, t) be an analytic vector 
field on the domain 

(z, t) e Vy X Vt 

Vy^izeC^-.lz- x„|oo < ry} 

Vt = {t e C : mt)\ < n} 

Then y(t) satisfying y' = g{y, t), j/(0) = x„ G M'^ is an analytic function of t on the 
domain 

V=\teC:\t\ <mm(-pl-,rt 

where \\g\\r :== sup(^_t)gx)„xX)t Iffl^'OU- 

We can now combine the estimates of Corollary 3 and Lemma 4 to determine a 
bound on p, and hence on the error in the numerically computed modified energy. 

Theorem 5 (Numerical modified energy). Let H(jp, q) he analytic in its arguments, 
and let Pn,qm Pn be computed by the algorithm (2). Then for each n there exists 
an m* such that we have the error bound 

|i:?(Pn,g„,t„)-T„",.^„,.| < ^expi~C2^J^] (kn|||p||p+|p„|||q||p + ||^||p), 

where C2 < 2.14707 (and p = 0.6835(5/||/||p;. 

The error bound in Theorem 5 shows that we are able to track the modified energy 
exponentially accurately. Moreover the bound displays the same dependency on the 
parameters h, fy and \\f\\v as Theorem 1. 

The bounding-constant C2 is however smaller than the 27r/e found in the proof 
of Theorem 1. It is unclear to us if this is due to the proof techniques applied, or if 
it is an actual weakness of the extrapolation method when applied to estimate the 
derivatives and thus the modified energy. 

6. Numerical experiments. We have implemented the extrapolation algorithm 
using the Arprec multiple-precision library in order to avoid pollution by round-off 
errors and to be able to verify the theory to high accuracy. ^ We set the precision to 
120 digits. Most experiments are done using the standard Stormer-Verlet scheme 
(also known as the leap frog scheme). All the experiments were also repeated with 
two fourth-order splitting schemes to check for dependence on splitting scheme 
coefficients. No noteworthy dependence was found, and we only present these results 
for the Kepler experiment. 

An early experiment verifying the exponentially small drift in modified energy 
was done by Benettin and Giorgilli [2] who used a Hamiltonian of the form H = 
5 {Pi +P2) + U{qf + g|) with the potential function U vanishing fast as its argument 
becomes large. In this case, the exponentially small effects can be observed directly 
because methods of the form (1) preserve the energy H exactly when U is identically 
zero. To carry out this experiment the initial values j/o = (pi(0)jP2(0),qi(0), (72(0)) 



^The Arprec library is available from http://crd.lbl.gov/~dhbailey/mpdist/. The CH — h 
source code for our experiments can be downloaded from http://wwwl.inaths.le6ds.ac.uk/ 
~j it se/sof tware.html 
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100 



Figure 1 . Drift in modified energy for the pendulum as a function 
of step size h for various values of m. The thick line indicates the 
limit. 



are then chosen so that U vanishes and that the trajectory passes close to {qi, q-z) = 
(0,0) before ending at some point ijt = {pi{T),p2{T),qi{T),q2{T)) where again U 
vanishes. Carrying out that simulation the difference |iJ(yo) ^ ^(z/t)| is observed 
to be 0{exp{—C/h)) where C is some unspecified constant. 

We repeated this experiment using our method, and found that in this case it had 
zero error so the experiments we consider will not have this type of Hamiltonian. 
This matter warrants further investigations, but we have not pursued these in this 
paper. 



6.1. The pendulum. In this experiment we apply the Stormcr-Vciiet method to 
the pendulum, which has Hamiltonian H = p^/2 — cos{q), integrated over the time 
interval [0,100]. Figure 1 reports the drift in the modified energy computed us- 
ing T„i_,„ for TO — 2,5,10,15,20,30,40. Here, and in the other plots, the drift is 
defined as the difference between the maximum of the modified energy over the 
integration interval and its minimum. The figure suggests that the approxima- 
tions Tm,m converge for this problem. The limit is indicated by the thicker line in 
the left part of the plot, which shows that the drift follows the exp(— c//i) behaviour 
predicted by the theory. 

The initial value for the experiment reported in Figure 1 is q(0) = and p(0) = 1. 
Next we study the effect of varying the initial condition. The result is shown in 
Figure 2. The Stormer-Verlet method shows improved energy preservation near 
the equilibrium point, revealing the exp(— c//i||/||x)) dependency on step size and 
on ll/llx' which decreases as p — > 0. 
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Figure 2. Drift in modified energy for the pendulum as a function 
of the initial momentum p(0). 



6.2. Kepler problem. The Kepler problem for one particle in a central force field 
is given by the Hamiltonian 

We integrate over the interval [0, 100] starting from the point 

Pi = 0, qi = 1 - ecc. 



P2 



1 



1 



92 



0. 



where < ecc < 1 is the eccentricity of the orbit. 

The left plot of Figure 3 shows the theoretical exp(— c//i) behavior. In contrast 
with the pendulum, for this problem the T„ „ do not converge as m — > oo, but the 
sequence has to be truncated at a suitably chosen point. To find the optimal m, we 
approximate the error in the ?7i-th estimate as 

We compute this estimate for rn = 2, 3, . . . , 200 and select the value of m for which 
the estimated error is minimized. This procedure recovers the expected exponential 
behaviour. 

The right plot of Figure 3 shows the dependence of the drift in the modified 
energy on the eccentricity of the orbit. Almost circular orbits with a low eccentric- 
ity show much better preservation of the energy than highly elliptical orbits. An 
instability occurs at a critical eccentricity which depends on the step size. This can 
be explained by the fact the the topology of the energy levels of the modified energy 
changes with /i, thus leading to unbounded trajectories and instability. 

We used the second-order Stormer-Verlet method to produce Figure 3. We ran 
the experiments again with two fourth-order splitting methods: Yoshida's scheme 
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Figure 3. Drift in modified energy for the Kepler problem as a 
fmiction of step size h for eccentricity ecc = 0.6 (left plot), and as 
a function of eccentricity (right plot). 



based on extrapolation [20] and a fourth-order scheme due to Blanes and Moan [3] . 
The last scheme is optimized for problems of the type we have considered. It has 
very small error coefficients at the cost of many stages, leading to coordinate errors 
which are typically three orders of magnitude smaller than Yoshida's method at 
the same computational cost. The plots for the drift in modified energy of both 
Yoshida's method and the Blanes-Moan method look the same as for the Stormer- 
Verlet method. In particular, the constant c in exp(— c//i) is the same. However, 
the drift in the modified energy for the Stormer-Verlet method is approximately a 
factor of three smaller than Yoshida's method and a factor of four smaller than the 
method of Blanes and Moan. 

The left plot in figure 4 illustrates for several different step sizes how the modified 
energy varies. There are peaks when the particle is near the singularity at the origin. 
The crucial point is that the energy essentially recovers its value after this point 
before another close encounter. The plot on the right compares the three different 
methods. It is seen that the methods give rather different results, even though the 
maximal variation in the modified energy is almost the same for the methods. The 
Blanes-Moan method seems to preserve the modified energy better after the close 
encounter, which might indicate a special advantage of this method when applied to 
the Kepler problem. If, however, the time steps are scaled so that the computational 
cost is the same for the three methods, the Stormer-Verlet method will preserve 
the modified energy better than the high-order methods. 

Figure 5 shows the accumulated change in energy, \H{po,qo,to) — H{pn, g„, t„)|, 
and the instantaneous change in energy, |7J(p„_i, g„_i, t„_i) — iJ(p„, g„, t„)|, to- 
gether with the optimal m found by the error estimate (6). The graph shows that 
near the singularity quite a high order m (which we bounded by 200) is used. This 
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Figure 4. The difference between modified energy at a given time 
and tlie initial modified energy for the Kepler method with eccen- 
tricity ecc = 0.6. The left plot shows the results for the Stormer- 
Verlet method for different step sizes, while the right plot shows 
the results for different methods. All methods are run with step 
size /i = O.f , so Stormer-Verlet does considerably less work. 



indicates that information from the smooth parts of the trajectory is used near 
singularities, and that it might be important to use very high order approximations 
to get a clear picture of the drift. The graph also indicates that the algorithm can 
track instantaneous changes in energy. 

Away from the parts of the orbit where the singularity at the origin is approached 
most closely, a lower value of m suffices. It is thus useful to find a more efficient 
method for finding the optimal m instead of computing the error estimate for all m 
up to some large value (here, 200). We found good results with the following ad-hoc 
termination criterion: compute the error estimate (6) for all m up to the first value 
of m for which 



max 

j— rn— 11,... ,m~ 



m 



3,0 



T, 



i-ij-i 



< 



IT, 



],3 



T 



J-lj-ll 



'1 j— m— 10,...,m 

and then choose the m with the minimal error estimate. The plots produced by 
this criterion are nearly indistinguishable from the plots produced when all m up 
to 200 are considered. 



6.3. Henon Heiles system 

by 

H ' 



The Hamiltonian of the Henon-Heiles system is given 



1 



ql 



Hq2 



U) 



Skeel et al. investigated the theoretical exp(— c//i) behaviour of the drift in the 
modified energy for this system, and report that the results are "less convincing" [5, 
§2.4]. We revisit this problem, using instead the extrapolation method to achieve 
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Figure 5. The change per step in the modified energy (dashed) 
and the accumlated change (sohd), and the optimal order m 
(capped by 200, marked by 'x'). The axis for the energies is on the 
left, which the right axis is for the order m. This is for the Stormer- 
Verlet method applied to the Kepler problem with h = 1/20 and 
ecc = 0.6. 




Figure 6. Drift in modified energy for the Stormer-Verlet method 
applied to the Henon-Heiles problem as a function of step size h 
with initial condition pi(0) = 0.1 (left plot), and as a function of pi 
for step sizes h — 0.25, 0.1, 0.05 (right plot). 
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arbitrary high orders. Figure 6 shows that the expected exponential smahness is 
indeed present, and that there is no problem in using the algorithm other than 
allowing for large values of to (we again capped m at 200). The effect of round-off 
error becomes visible when 1/h exceeds 20; remember that all computations are 
done with 120 digits. 

The right plot shows how the maximal deviation of the modified energy changes 
as the initial condition for pi is varied; the initial conditions for the other variables 
are fixed as qi = q2 = P2 = 0. This plot shows that there is no abrupt change in 
energy preservation when moving from regular, integrablc motions (the region with 
energy H < 1/12 or, equivalcntly, pi < 1/a/6 « 0.4) to the chaotic regime of phase 
space (where H > 1/12). 

We also applied the fourth-order methods due to Yoshida and Blanes and Moan 
to this problem, with the same results as for the Kepler problems: the Stormer- 
Verlet method shows slightly better energy preservation, but the value of c in the 
exp(— c//i) dependence is the same. 

7. Conclusions. We have supplied rigorous estimates for a numerical algorithm 
that computes the modified energy for methods based on operator splitting of Hamil- 
tonian systems. The estimate shows that the procedure can recover exponentially 
small estimates, known to exist theoretically. The estimates exhibit the same depen- 
dence on the important parameters fy, h and ll/Hx), and can therefore in principle 
be used to extract their values from simulations. When comparing different split- 
ting algorithms, it seems that in the limit /i -> the exponential remainder term 
only weakly depends on the method coefficients. Thus when considering the addi- 
tional cost of optimized, many-stage, methods these will have a larger drift than 
the second-order Stormer-Verlet algorithm. In other words, when it comes to pre- 
serving the modified energy, cheap, low-order methods are preferable. Although 
we have not considered ODEs originating from Hamiltonian semidiscretization of 
PDEs it seems likely that for long time simulations a low-order method such as 
Stormer-Verlet is preferable if energy preservation is important. 

Appendix. 

Proof of Lemma 2. Without loss of generality we assume that n — Q. 
By representing (4) as a contour integral we have 

1 ^ r {-iy+^{m\f f 1 i_\ 

2Tri^Jhj{m-jy.{m + jy.\z-jh z + jh j 



J- 



=i-J'< 



where the contour 7 includes the points —mh, . . . ,mh and excludes singularities 
of y, as sketched in Figure 2. 

The derivative is given by y'{Q) — 5;^ / ^ J2 ^ , so the error in the approximation 
becomes 

E„,iym = A„y(0) - 2/'(0) = -^ / K„,iz)y{z)dz, (7) 



2ni j^ 



where the kernel is defined by 

(^1)™+1(to!)2/i2" 



Km{z) 



z'^{z'^ ~ h^){z^ ~ {2hY) • • • (z2 - (m/i)2) ■ 
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-s{z) 



X 



X 




Xp -(m-l)ft 

---o e e 

^"'^ -(m - 2)h 



-o---p--k>- 

-h ' h 




X 



X 



Figure 7. The contour of integration used in the proof of Lemma 2. 



Along the curve 7, the kernel Km achieves its maximum in modulus &t z — ip, and 
the maximum is 



\Kmiip)\ = 



(m!)2/i2" 



1 



p2(p2 + h^)... (p2 + (^/j)2) ^2(1 + £i) . . . (1 + pL^) 



^ TT I P 

= phsinHY)J}^X^^^ 

where the last equality follows from n?li(-'^ + (h)^ ^ ^ ^ sinh(7rp/ft,). The product 
can be bounded as 



log Y[ (1 

j—m-\-l ^ 



(hjy 



E i°g 1 



j—m-\-l 

< log ( 1 



P 



{hjr 

r,2 



(hx 



dx < 



h^m^ 



yielding UZm+ii^ + TKw) < exp(j^), and thus 



\Kmiip)\ < , „■ , /TTpx exp 



Since the length of the contour is 2'Kp+Amh, the error expression (7) can be bounded 



as 



\Dmy{0)-y'm< 



(7rp + 2m/i)exp(^^j 



p/isinh(^) "^"^ 

where ||y||p := supicj^^-ji^p |?;(z)|oo- This upper bound is minimized by choosing m 



so that 2^(7rp + 2m/i) exp{p^/h^m,) vanishes, i.e. m ^ j^. This gives the bound 



iD^yio) - .'(0)1 < ^(:l?f/,f ||.t < cZ-^^'^^f^M, 



h sinh I 



/i2 



for some constant Ci > 0. 



n 



14 



PER CHRISTIAN MOAN AND JITSE NIESEN 



Proof of Corollary 3. This follows from 

\H{p.n,qn,tn)-T^n',,n'\ 

< k\qn\l\\p' - Dra'VWp + jKllllg' - An'^llp + \\\^ ~ Dr,-.'Mp 



< 



Cipexp(-f) 



(knlillpllp + Klillgllp + ll/^llp) 



where the last inequality follows from Lemma 2. 



a 



Proof of Lemma 4- We prove this by Picard iteration: set xi — x„ and iterate 
Xk+i{t) = Xn + Jg dixkis), s) ds. Fix t G 2?t, and assume at first that rt is sufficiently 
large. For Xk+i,Xk G Vy 

\g{xk+i,t) - g{xk,t)\oc 

—g(xk+i + s{xk - Xk+i),t] ds 
f^ f g{xk+i + zixk - Xk+i), t) 

Jo /|2-s|=_R 
fl 




f 

2^ 



dz ds 



1 

2^ 



g{xk+i + s{xk - Xk+i) + w(xk - Xk+i),t) 



dw ds 



10 J\w\=R <" 

The radius R is restricted by the requirement that the argument of g lies within Dy 



\xk + l + s{xk - Xk+l) - Xn + w(xk - Xfe+i)|oo 

< \xk+l + s{xk - Xk+l) - a^nloo + r\{xk - Xk + l))\oo 

< \t\\\g\\r + R\xk -Xk+i\oo < ry 
by using 

\xk+l + s{xk - Xk+l) - Xn)\oo 

(1 - s) \g{xk{T), r)|oo + s \g{xk-\{T), r)|oo dr < \t\\\g\\, 



< sup 

\r\ = \t\ 



We may therefore choose 



R = T]] 



ItMi 



\Xk+l Xk\oo 

which gives the supremum-norm Lipschitz bound 

ll.9l|r 



-, 0<ri<l 



\g{xk+i,t) - g{xk,t)\^ < 



\'^k-\-l •^k\oo- 



Vi-Ty- \t\\\9\\r) 

Let Afc+i(|i|) = supj^i 1(1 |a;fc+i(T) — a;fe(r)|oo, then the Picard iteration xi — Xn, 
Xk+l = Xn + /q g{xk{s), s) ds, converges if A^ — > as fc — > oo, with 



Ak+iit) < f 
Jo 



V{ry-s\\9\\r) 



Ak{s)ds, Ai(i) = |i|||5|| 
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Introducing the generating function G(ji, \t\) — X]fc>i m'^^A; we have 

G{i,,t)<fi\t\\\g\\r + fi f , "^";:„„ G{^i,s)ds. 
Jo Vvy — s\\J\\p) 

Since the terms in this inequahty are positive, an upper bound is the solution of 

^^^^ = ^^\\9\\r + , "^I'liii II , g+(a^, \t\), g+(m, \t\ = 0) = 0, 

d\t\ Viry - \t\\\9\\r) 



G+M\) ^J^ffl^mk)-'^' (,_\tML' 

Because G"*" is analytic in ^ around /i = 1 provided |i| < jzhi the sequence Afe(|i|) 
converges uniformly to zero and hence Xk(t) converges uniformly to the solution. 
Since each iterate Xk+i{t) ~ a;„ + J^ f{xk{s),s)ds is analytic in t G {t G C : 
\t\ < niinj^T-, ?'(}} the uniform convergence gives by Weierstrass theorem that 
y{t) — Xao(t) is analytic in this domain as well. D 

Proof of Theorem 5. In Theorem 1 we take g = f, thus \\g\\r < t^II/IIp with / 
analytic in V . This gives that the y{t) is analytic in the domain 



We find that the bound is optimized by picking fy — rjS where rj — '^ ^— . Thus 
we may take p ~ -rw— < 0-6835 n A in Corollary 3 giving the exponentially small 



bound (t = nh) 

\H{pn,qn,t) -r,„.,™.| < ^exp f-G2 ^,, ,, j (k«|i||p||p + bn|i||g||p + Pllp), 
where G2 = 0.68347r. D 
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